﻿using MathNet.Numerics;
using SeeSharpTools.JY.DSP.Fundamental;
using SeeSharpTools.JY.DSP.Utility.Filter1D;
using Sunny.UI;
using System;
using System.Drawing;
using System.Linq;

namespace SeeSharpTools.Demo
{
    public partial class FormFIRDelay : UIPage
    {
        public FormFIRDelay()
        {
            InitializeComponent();
        }

        private void FormFIRDelay_Load(object sender, EventArgs e)
        {
            // 生成含噪声的心电图信号
            double fs = 500;
            int length = 500;
            double[] x = ECG(length);
            double[] noise = new double[length];
            Generation.UniformWhiteNoise(ref noise, 0.25);
            x = x.Zip(noise, (a1, b1) => a1 + b1).ToArray();

            // 设计截止频率为 75 Hz 的 70 阶低通 FIR 加窗滤波器，并验证线性相位滤波器的群延迟为常量
            double[] wn = new double[] { 70 / (fs / 2) };
            double[] b = FIRDesign.Window(70, wn, bandType: FIRBandType.Lowpass);
            var (gd, f) = FIRAnalysis.GroupDelay(b, fs);
            easyChartX_GrpDelay.Plot(f, gd);
            easyChartX_GrpDelay.Series[0].Name = "群延迟";
            easyChartX_GrpDelay.AxisX.Title = "频率(Hz)";
            easyChartX_GrpDelay.AxisY.Title = "延迟(Samples)";

            // 不进行延迟补偿，直接对原始信号进行滤波
            double[] y1 = FIRFiltering.FIRFilter(b, x);

            // 在原始信号末尾追加延迟个数的零，然后进行滤波，
            // 移除滤波后信号最前端延迟个数的样本点，保证输出信号和输入信号长度相等，消除滤波器瞬变
            int delay = (int)gd.Average();
            double[] y2 = FIRFiltering.FIRFilter(b, x.Concat(new double[delay]).ToArray());
            y2 = y2.Skip(delay).ToArray();

            double[,] data = new double[3, x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                data[0, i] = x[i];
                data[1, i] = y1[i];
                data[2, i] = y2[i];
            }
            easyChartX_Filtering.Plot(data, 0, 1 / fs);
            easyChartX_Filtering.Series[0].Name = "原始信号";
            easyChartX_Filtering.Series[0].Color = Color.Orange;
            easyChartX_Filtering.Series[1].Name = "不补偿延迟滤波";
            easyChartX_Filtering.Series[2].Name = "补偿延迟滤波";
            easyChartX_Filtering.AxisX.Title = "时间";
            easyChartX_Filtering.AxisY.Title = "幅度";
        }

        /// <summary>
        /// 生成指定长度的分段线性的ECG信号
        /// </summary>
        /// <param name="length">信号长度</param>
        /// <returns>输出的ECG信号</returns>
        public static double[] ECG(int length)
        {
            double[] a0 = new double[] { 0, 1, 40, 1, 0, -34, 118, -99, 0, 2, 21, 2, 0, 0, 0 }; // Template
            double[] d0 = new double[] { 0, 27, 59, 91, 131, 141, 163, 185, 195, 275, 307, 339, 357, 390, 440 };
            double[] a = a0.Select(num => num / a0.Max()).ToArray();
            int[] d = new int[d0.Length];
            for (int i = 0; i < d0.Length; i++)
            {
                d[i] = (int)Math.Round(d0[i] * length / d0[14]); // Scale them to fit in length
            }
            d[14] = length;

            double[] x = new double[length];
            int[] m;
            for (int i = 0; i < 14; i++)
            {
                m = Generate.LinearRangeInt32(d[i], 1, d[i + 1] - 1);
                double slope = (a[i + 1] - a[i]) / (d[i + 1] - d[i]);
                for (int j = 0; j < m.Length; j++)
                {
                    x[m[j]] = a[i] + slope * (m[j] - d[i]);
                }
            }
            return x;
        }
    }
}